Consignes

Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :

La bioinfo c'est : <code>MERVEILLEUX</code>.

N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.

Introduction

Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Complete Genome Sequences of 13 Bacillus subtilis Soil Isolates for Studying Secondary Metabolite Diversity

Analyses

Organisation de votre espace de travail

La première étape est la création de differents répertoires de travail:

FASTQ : pour stocker les fichiers téléchargés
QC : Pour stocker les résultats des analyses qualités des fichiers fastq avant et après nettoyage
CLEANING : Pour stocker les résultats du nettoyage des fichiers fastq
MAPPING : pour le génome de référence, son annotation et les résultats des alignements

# Dans le répertoire "EvaluationM4M5-main" on crée les différents répertoires de travail
mkdir -p ~/shared/home/hbelcram/GITRepository/Module_4_5_Outils_Haut_debit/Module_5_transcript/EvaluationM4M5-main/CLEANING
mkdir -p ~/shared/home/hbelcram/GITRepository/Module_4_5_Outils_Haut_debit/Module_5_transcript/EvaluationM4M5-main/FASTQ
mkdir -p ~/shared/home/hbelcram/GITRepository/Module_4_5_Outils_Haut_debit/Module_5_transcript/EvaluationM4M5-main/MAPPING
mkdir -p ~/shared/home/hbelcram/GITRepository/Module_4_5_Outils_Haut_debit/Module_5_transcript/EvaluationM4M5-main/QC

Téléchargement des données brutes

Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]

# Chargement du module sra-toos pour utiliser la fonction fasterq-dump de téléchargement des fichiers fastq
module load  sra-tools
# Information sur la version
fasterq-dump --version
#"fasterq-dump" version 2.10.3
# Lancement du téléchargement:
srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir FASTQ

Combien de reads sont présents dans les fichiers R1 et R2 ?

le téléchagement affiche un total de 14 132 110 reads. Mais je procède à la compression des fichiers fastq pour compter les reads.

# compression des fichiers fastq en .gz
# Dans le répertoire FASTQ je lance la commande suivante:
srun gzip *.fastq
# Pour vérifier que la taille des fichiers
ls -lh 
# Comptage du nombre de ligne par fichier fastq.gz:
# Pour le fichier 1:
echo $(zcat /shared/home/hbelcram/GITRepository/Module_4_5_Outils_Haut_debit/Module_5_transcript/EvaluationM4M5-main/FASTQ/SRR10390685_1.fastq.gz |wc -l)
# Pour le fichier 2:
echo $(zcat /shared/home/hbelcram/GITRepository/Module_4_5_Outils_Haut_debit/Module_5_transcript/EvaluationM4M5-main/FASTQ/SRR10390685_2.fastq.gz |wc -l)

Les fichiers FASTQ contiennent 7066055 reads chacun:

28264220/4 soit 7066055 reads pour SRR10390685_1.fastq.gz
28264220/4 soit 7066055 reads pour SRR10390685_2.fastq.gz

Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

# Téléchargement du génome de référence dans le répertoire MAPPING:

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz

Quelle est la taille de ce génome ?

# Décompression du fichier génome
srun gunzip GCF_000009045.1_ASM904v1_genomic.fna.gz
# Pour la taille on enlève la ligne du nom (grep -v) et on compte le nombre de caractères (wc -m):
grep -v '>' /shared/home/hbelcram/GITRepository/Module_4_5_Outils_Haut_debit/Module_5_transcript/EvaluationM4M5-main/Data/GCF_000009045.1_ASM904v1_genomic.fna | wc -m

La taille de ce génome est de 4 268 302 paires de bases (pb) soit environ 4,27Mb. Sur le site du NCBI on trouve l’information d’une taille de 4,215,606 pb (https://www.ncbi.nlm.nih.gov/assembly/GCF_000009045.1/)

Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

# Télechargement du fichier génome de Bacillus dans le répertoire MAPPING:
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz

# décompression du fichier d'annotation du génome:
srun gunzip GCF_000009045.1_ASM904v1_genomic.gff.gz

Combien de gènes sont connus pour ce génome ? Le nonbre de gènes connus pour ce génome est de 4,536 gènes suivant le site NCBI (https://www.ncbi.nlm.nih.gov/genome/?term=txid1423[orgn])

# Comptage du nombre de gènes dans le ficher d'annotation:
grep 'ID=gene' GCF_000009045.1_ASM904v1_genomic.gff | wc -l

4536 gènes sont recensés dans le fichier d’annotation, ce qui correspond à l’information trouvé sur le ste NCBI.

Contrôle qualité

Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit

# chargement de l'outil d'analyse qualité: fastqc dans le répertoire EvaluationM4M5-main et dépôt des résultats dans le répertoire QC
module load fastqc
fastqc --version
#FastQC v0.11.9
# lancement de l'outil fastqc sur les fichers fastq
srun --cpus-per-task 4 fastqc ../FASTQ/SRR10390685_1.fastq.gz -o ../QC/ -t 4
srun --cpus-per-task 4 fastqc ../FASTQ/SRR10390685_2.fastq.gz -o ../QC/ -t 4

La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

[X ] Oui
[] Non

les reads ont des notes de qualités supérieures à 28 quand on regarde les résultats obtenus avec l’outil fastqc. Ce qui veux dire une erreur de lecture tout les 1000 bases. On observe dans le rapport pour les reads du fichiers SRR10390685_2, un problème d’éllimination des adaptateurs sur les reads les plus longs (Adapter content). Le “Quality” score ne dépasse pas 36.

Lien vers les rapports FastQC

Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

[] Oui
[x] Non

La distribution montre une majorité des reads autour de 150pb. Ceci pourrait provenir de la taille des fragments lors de la construction des librairies ou d’un probléme autre. les tailles ne sont pas celles visées dans l’article autour de 250pb.

Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?

La profondeur de séquençage : nombre de reads x tailles des reads / taille du génome soit: 7066055 x 150 / 4268302 soit 248,32 X. La profondeur de séquençage est d’environ 250X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.

Les paramètres suivants ont été choisis : Utilisation de l’outil fastp avec les paramètres principaux suivant: - Sélection des reads avec une note minimale de 30 (cut_mean_quality 30) - La taille minimale des reads à 100pb (length_required 100)

# Chargement de l'outil fastp pour filter les fichiers fastq:
module load fastp
fastp --version
#fastp 0.20.0

# Lancement  de fastp pour le néttoyage:
srun --cpus-per-task 4 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 4 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json CLEANING/fastp.json

Ces paramètres ont permis de conserver 13 554 096 reads pairés, soit une perte de 4,09% % des reads bruts.

Vérification des reads nettoyés avec fastp par l’outil fastqc

# Utilisation de fastqc sur les fichiers néttoyés:
# chargement de l'outil d'analyse qualité: fastqc
module load fastqc
fastqc --version
#FastQC v0.11.9
# lancement de l'outil fastqc sur les fichers fastq
srun --cpus-per-task 4 fastqc ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz -o ../QC/ -t 4
srun --cpus-per-task 4 fastqc ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz -o ../QC/ -t 4

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].

# Dans le répertoire MAPPING
# Chargement de l'outil de mapping bwa
module load bwa
# indexation du fichier séquence du génome
bwa index GCF_000009045.1_ASM904v1_genomic.fna
# bwa version?

Lancement du mapping des reads avec bwa:

# Alignement des reads avec bwa
srun --cpus-per-task=32 bwa mem GCF_000009045.1_ASM904v1_genomic.fna ../CLEANING/SRR8082143_1.cleaned_filtered.fastq.gz ../CLEANING/SRR8082143_2.cleaned_filtered.fastq.gz -t 32 > SRR8082143_on_ASM904v1_genomic.sam

Transformation du fichier .sam en fichier .bam à l’aide de samtools pour la suite des analyses:

# chargement de l'outil samtools
module load samtools
# recherche de la version
samtools --version
# samtools 1.10
# Using htslib 1.10.2
# Utilisation de samtools:
srun --cpus-per-task=8 samtools view --threads 8 SRR8082143_on_ASM904v1_genomic.sam -b > SRR8082143_on_ASM904v1_genomic.bam

Combien de reads ne sont pas mappés ?

# Quelques statistiques sur l'alignement des reads avec samtools idxstats et samtools flagstat:
srun samtools idxstats SRR10390685_on_ASM904v1.sort.bam > SRR10390685_on_ASM904v1.sort.bam.idxstats
srun samtools flagstat SRR10390685_on_ASM904v1.sort.bam > SRR10390685_on_ASM904v1.sort.bam.flagstat

# Visualisation des sorties idxstats et flagstats
cat SRR10390685_on_ASM904v1.sort.bam.flagstat
cat SRR10390685_on_ASM904v1.sort.bam.idxstats

Résultats:

2252171 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 2905 + 0 supplementary 0 + 0 duplicates 2129814 + 0 mapped (94.57% : N/A) 2249266 + 0 paired in sequencing 1124633 + 0 read1 1124633 + 0 read2 2116428 + 0 properly paired (94.09% : N/A) 2120300 + 0 with itself and mate mapped 6609 + 0 singletons (0.29% : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5)

NC_000964.3 4215606 2129814 6609 * 0 0 115748

6609 reads ne sont pas mappés.

Croisement de données

Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:

# Recherche des reads couvrant le gene trmFNF avec l'outil bedtools:
module load bedtools
srun bedtools --version
#bedtools v2.29.2
# Conversion du fichier .sam en .bed: SRR10390685_on_ASM904v1.sort.bed pour l'outil:
srun bedtools bamtobed -i SRR10390685_on_ASM904v1.sort.bam > SRR10390685_on_ASM904v1.sort.bed

# Recherche des reads couvrants le gène trmNF à 50%:
srun bedtools intersect -a SRR10390685_on_ASM904v1.sort.bed -b trmNF_gene.bed -f 0.50 > SRR10390685_on_trmNF_gene.bed
# Comptage des lignes (équivalent au nombre de reads) dans le fichier de sortie:
wc -l SRR10390685_on_ASM904v1.sort.bed
# 469

469 reads chevauchent le gène d’intérêt.

Visualisation

Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.

# Recherche des reads couvrant le gene trmFNF:
module load bedtools
bedtools --version
#bedtools v2.29.2
# Convertion du fichier BAM au format BED pour utiliser bedtools.
srun bedtools bamtobed -i SRR10390685_on_ASM904v1.sort.bam > SRR10390685_on_ASM904v1.sort.bed

# Recherche des reads couvrant le gene trmNF au moins de 50%:.
srun bedtools intersect -a SRR10390685_on_ASM904v1.sort.bed -b trmNF_gene.bed -f 0.50 > SRR10390685_on_trmNF_gene.bed

#![Visualisation des reads sur le gène trmNF]()

#![Zoom des reads sur le gène trmNF]()

Visualisation des reads sur le gène trmNF avec le visualisateur ‘genomeview’ En entrée:

le génome
l’annotation du génome
les 469 reads trouvés pour le gène

On obtient:

Visualisation des reads sur le gène trmNF

Zoom des reads sur le gène trmNF

References

1. toolkit NS. NCBI sra toolkit. NCBI, GitHub repository. 2019.

2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

3. Zhou Y, Chen Y, Chen S, Gu J. Fastp: An ultra-fast all-in-one fastq preprocessor. Bioinformatics. 2018;34:i884–90. doi:10.1093/bioinformatics/bty560.

4. Li H. Aligning sequence reads, clone sequences and assembly contigs with bwa-mem. arXiv preprint arXiv:13033997. 2013.

5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and samtools. Bioinformatics. 2009;25:2078–9.

6. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

7. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (igv): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France

Examen final Modules 4 et 5

DUBii 2021

Harry Belcram

14 April, 2021